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Abstract 

We address a detailed non-perturbative numerical study of the scalar theory on the fuzzy 
sphere. We use a novel algorithm which strongly reduces the correlation problems in the 
matrix update process, and allows the investigation of different regimes of the model in a 
precise and reliable way. We study the modes associated to different momenta and the role 
they play in the "striped phase" , pointing out a consistent interpretation which is corroborated 
by our data, and which sheds further light on the results obtained in some previous works. 
Next, we test a quantitative, non-trivial theoretical prediction for this model, which has been 
formulated in the literature: The existence of an eigenvalue sector characterised by a precise 
probability density, and the emergence of the phase transition associated with the opening of a 
gap around the origin in the eigenvalue distribution. The theoretical predictions are confirmed 
by our numerical results. Finally, we propose a possible method to detect numerically the 
non-commutative anomaly predicted in a one-loop perturbative analysis of the model, which 
is expected to induce a distortion of the dispersion relation on the fuzzy sphere. 



1 Introduction 



The study of quantum field theory in non-commutative spaces has attracted considerable attention 
over the last years [1-4]. This research area has a long history, since the possibility of a quantised 
structure of spacetime at short distances was first mentioned as early as in the 1930's in some 
correspondence among Heisenberg, Peierls, Pauli and Oppenheimer [5, 6], and in the papers 
published by Snyder [7], by Yang [8] and by Moyal [9] during the 1940's. Although the original 
motivation to use non-commutativity as a tool to regularise QFT was soon frustrated — while the 
renormalisation approach proved to be a successful method to handle the divergences encountered 
in the formulation in commutative spacetime — , non-commutative spaces have attracted renewed 
interest in more recent years, with the application of this formalism to solid-state physics, to the 
problem of the quantum Hall effect [10], to the study of MHD waves in astrophysics [11], and with 
the discovery of the relevance of such spaces to string theory [12-18] and to a possible quantum 
theory of gravity [19,20]. 

Groenewold-Moyal R^ spaces are among the most extensively studied non-commutative spaces. 
The properties of QFT defined in these spaces — including those related to renormalisabil- 
ity, causality, non-locality, Poincare invariance et ccetera — have been addressed in several 
works [21-44]; a lattice-like regularised formulation has also allowed to investigate numerically 
various aspects of these models [45-54]. In particular, one of the most interesting — albeit trou- 
blesome — features is the fact that the effective action describing QFT in a Groenewold-Moyal 
space is divergent when the external momentum along the non-commutative directions vanishes: 
this effect arises from the integration of the high-energy modes in non-planar loop diagrams, and 
is henceforth called "ultra- violet /infra-red (UV/IR) mixing". 

Another class of non-commutative spaces is given by fuzzy spaces: they are built approxi- 
mating the infinite-dimensional algebra of functions on some particular manifold by means of a 
finite-dimensional algebra of matrices. Under some conditions, this construction is possible for 
even-dimensional co-adjoint orbits of Lie groups which are symplectic manifolds — see [55-68] 
and references therein. In particular, co-adjoint orbits of semi-simple Lie groups are adjoint or- 
bits; examples include the CP n complex projective spaces. The most-widely known example of 
a fuzzy space is the fuzzy two-sphere Sp [55], built truncating the algebra of functions on the 
commutative sphere S 2 to a maximum angular momentum / max . The fuzzy sphere depends on 
two parameters: the matrix size ./V = / max + 1 and the radius R; the commutative sphere and the 
non-commutative plane can be obtained in different limits of N and R. 

A one-loop perturbative calculation shows that, for every finite N, QFT on the fuzzy sphere 
is finite; furthermore, the theory is not affected by the UV/IR mixing problem [69] (although — 
due to a non-commutative anomaly — the latter re-emerges once a double-scaling limit is taken, 
in which the fuzzy sphere goes over to the non-commutative plane Rg). This feature, as well as 
the fact that the fuzzy approach explicitly preserves the symmetries of the original manifold for 
any value of N and allows a well-defined treatment of the topological properties [70-84] , has led 
to suggest the fuzzy space as a potentially interesting candidate for regularisation of quantum 
field theory. 

As a matter of fact, QFT on the fuzzy sphere is mathematically well-defined and finite [85], 
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and the formulation is amenable to a non-perturbative approach and to numerical studies using 
Monte Carlo simulations, with fields represented as finite-dimensional matrices. This approach 
has been followed in various recent works [86-92]. 

In the present paper we address a detailed Monte Carlo study of the <I> 4 scalar field theory 
on the fuzzy sphere; among other issues of interest, this simple model provides a laboratory to 
test the possibility to use the fuzzy space approach as a potential regularisation scheme for more 
realistic field theories. As it concerns the practical implementation of Monte Carlo simulations of 
the model, we shall present a novel algorithm, which reduces the autocorrelation time, combining 
overrelaxation steps with ergodic configuration updates. 

Preliminary results of this study have been presented in [93]. 

This manuscript has the following structure: in section [2] the theoretical framework under- 
lying the model is recalled, and the basic notations are introduced. In section [3] we discuss the 
implementation of the numerical simulations of the model, and present the results obtained with 
our algorithm. In section H] we comment on the implications of these results and on possible 
research perspectives. A technical discussion of the algorithm is presented in the Appendix [Al 

2 Review of the construction of the model 

A general discussion of the mathematical construction of fuzzy spaces can be found in many 
excellent articles and books, like, for instance [67]; for the scalar field theory on the fuzzy sphere 
Sp, we refer the reader to the detailed presentation in [85]. 

The basic idea is to replace the infinite-dimensional, commutative algebra of polynomials 
generated by the {xi}i=i,2,3 coordinates on the two-dimensional sphere XiXi = R 2 embedded 
in M 3 with a non-commutative algebra generated by {£i }i=i,2,3 operators satisfying a (trivially 
rescaled) su(2) Lie algebra. The latter can be realised using the Wigner-Jordan construction of 
the su(2) generators, restricting to the finite-dimensional (./V-dimensional) subspace of the Fock 
space generated by a pair of mutually commuting creation operators. Accordingly, the algebra 
of functions on the commutative sphere S 2 is replaced by the Matjv algebra, whose elements can 
be expanded into irreducible representations of su{2). 

The fuzzy sphere admits the commutative sphere S 2 and the Groenewold-Moyal plane M. 2 , as 
two different limits: the former is recovered for N — > oo with R fixed, whereas the latter can be 
obtained — at least locally — via a stereographic projection from a fixed point, in the double 
limit: N — > oo, R — ► oo, keeping R 2 /N fixed. 

The action for a massive, neutral, scalar field with quartic interactions on the fuzzy can be 
defined — according to the conventions used in [86] — as: 





where $ E Mat at is hermitian and can be expanded in the {1^ )m } polarisation tensor basis as: 
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The model can be quantised in the path integral approach, defining the expectation values of 
generic observables O = 0(ci ;7n ), which can be evaluated perturbatively, or estimated numerically 
from Monte Carlo simulations. 

The perturbative treatment of the theory is formulated via a proper definition of the Feynman 
rules — in particular, the interaction vertices are modified in a non-trivial way, which depends 
on N, and is consistent with the fusion rules. A careful one-loop analysis of this model [69] shows 
that the UV/IR mixing phenomenon does not occur on S F ; however, a non-commutative anomaly 
shows up, as a finite difference between planar and non-planar tadpole diagrams. 

This anomaly is expressed by a rotationally invariant, non-local contribution to the quantum 
effective action; it distorts the dispersion relation, and survives the limit to the commutative 
sphere. This seriously threatens the possibility to consider the fuzzy approach as a bona fide 
regularisation scheme for theories defined in a commutative space; however, according to [94], the 
problem may be overcome, redefining the interaction term in the matrix action with a normal- 
ordering prescription, which allows to cancel the undesired momentum-dependent quadratic terms 
in the effective action. An alternative, and more general, possibility would be to include rota- 
tionally symmetric higher derivative terms in the action, as suggested in [95]. 

On the other hand, when the Groenewold-Moyal plane limit is taken, the non-commutative 
anomaly reproduces the logarithmic divergence which is characteristic of the UV/IR mixing [28]. 

Fuzzy spaces were studied in [33] as a mean to regularise scalar field theory in non-commutative 
Rg spaces. There, it was shown that different phases can be distinguished in the large N limit, 
according to the form of the distribution of eigenvalues associated to the matrix <£. The argument 
goes as follows: For the free case, in the large N limit and assuming that the cut-off is much 
larger than the non-commutative scale, the leading contribution to the expectation values of even 
powers (2k) of the field come from planar diagrams. Next, one observes that, in the large N limit, 
a clustering property for the expectation values of products of integrals of the field holds, and 
implies that the measure is strongly localised. This allows to identify the eigenvalue distribution 
as the Wigner semi-circle law, which holds for gaussian random matrix models [96,97]. 

Next, one can generalise to include interactions: for a scalar theory with bare square mass 
m 2 and quartic interactions with coupling | in non-commutative space the resulting eigenvalue 
distribution is known exactly [98]: for m? larger than a critical value m 2 lit the eigenvalue density 
Pg(tp) has a connected support and is given by: 



where g' and {m') 2 are related to g and m 2 — see [33] for the details. On the contrary, for 
m 2 < m 2 iit , the eigenvalue distribution exhibits two disconnected peaks of finite width. 

Although the arguments underlying this derivation are expected to hold only approximately 
in the two-dimensional case (because the dominance of the planar diagrams over the non-planar 
ones is weaker than in four dimensions), one can check if the numerical results are consistent with 
the predicted features. 

The real scalar model with quartic interactions on the fuzzy sphere was studied numerically 
in [86,87], where it was shown that the model exhibits three different phases: a disordered 
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phase, in which the field typically fluctuates around zero; a uniform order phase, characterised by 
fluctuations around the broken-symmetry minima of the potential, and a non-uniform order phase, 
which was described as new, intermediate, phase, intrinsically related to the matrix nature of the 
fuzzy regularisation. The appearance of the latter was interpreted assuming that, in a certain 
parameter range, the kinetic contribution to the action might be negligible, and the dynamics of 
the system were effectively reduced to the framework of a pure potential model [99,100]. This 
phase was also described as analogous to the striped phase predicted in [30] for non-commutative 
Groenewold-Moyal spaces, and observed numerically in [51,52]. 

3 Numerical simulations 

The numerical approach to the model is completely straightforward, and, under many respects, 
analogous to the more conventional lattice settingl^J Expectation values of the observables are 
estimated from averages over finite ensembles of matrices {$}, characterised by a statistical 
weight which depends on the model dynamics; the algorithm generating the matrix ensemble is 
built combining overrelaxation steps [101,102] and canonical updates. 

For the configuration-updating process, different types of pseudo-random number generators 
were compared; the G05CAF generator of the NAG library has eventually been used in the 
production runs. 

For each choice of the (N, r, A, R) parameters, the autocorrelation time between elements in the 
thermalised matrix ensemble has been calculated using the auto-windowing procedure [108]; the 
expectation values of the various observables have been evaluated from ensembles of statistically 
uncorrelated matrices. The data analysis was done using standard techniques, and errorbars have 
been estimated using the jackknife method — see, for instance, [109]. 

In the r < regime, it is particularly interesting to study the behaviour associated with the 
various /-modes. The classical minima of the potential correspond to uniform distributions and 
obviously their physical content is purely described by the scalar (1=0) channel. Since we are 
dealing with a quantum model and the system size is finite, the ground state is actually unique, 
as quantum fluctuations allow finite-action tunneling events between the two minima. When r is 
negative in sign and large in modulus, the profile of the potential is very steep, and the typical 
matrix configurations lie in a close neighbourhood around the classical minima; the trace of $ 
allows to identify around which of the two minima the matrix is lying at a given Monte Carlo 
time. The "trace susceptibility" , defined as: 

X = <(trd>) 2 > - (|trd>|) 2 , (4) 

1 lt should be noted, however, that the fuzzy space approach to a quantum theory differs with respect to the 
lattice formulation in a number of aspects, which also have practical implications for the computer simulations. In 
particular, in lattice field theory there exist many efficient update methods (including, for instance, those described 
in [101-107]) which are based on the locality of the discretised action; these methods allow to strongly damp the 
autocorrelation among subsequent configurations in the Markov chain. Furthermore, parallel computation can 
often be implemented in a straightforward way. On the contrary, in the fuzzy setting, the dynamics of each degree 
of freedom is non-trivially entangled with each other's, and — in general — the implementation of analogous 
techniques is not trivial. 
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encodes the physical information about the fluctuations of tr$. When r is increased to values 
closer to zero, the trace susceptibility exhibits a peak, corresponding to a maximum in the 
quantum fluctuations. As usual, the location of a maximum in the susceptibility approximatel)H 
identifies the critical value where the "phase transition" to the disordered phase would occur, for 
an infinite system. 
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Figure 1: Behaviour of the trace susceptibility, as defined in eq. 

The model already exhibits good scaling properties for small matrix sizes N < 10. As an 
example, figure [T] shows the results obtained for the susceptibility, along the A = ^ line: the 
peak gets very pronounced as N is increased, and its location is not affected by strong finite- iV 
effects. 

Taking a closer look at the matrix ensembles and at their physical content in non-zero mo- 
menta allows to detect the non- uniform modes, and the role they play in the regime under 
consideration. A simple variational analysis shows that also the commutative sphere can admit 
non-uniform configurations characterised by a finite, negative total amount of actionH Numer- 
ically, the relevance of these configurations can be quantified through the average values of the 
square moduli of the q m coefficients for I > 1. The latter are indeed found to be non- vanishing 
for all of the N, A and r values which we investigated; again, in agreement with [86,87], we 
observed a parameter range where the expectation value of the I = 1 mode is larger than the 
scalar component. 

The transition from the uniform- to the non-uniform-order phase can be interpreted as an 
2 Modulo corrections due to the finiteness of iV and/or R. 

3 It is easy to verify analytically that — at least in some parameter ranges — even an axially symmetric (but 
non-uniform) configuration described in terms of the first Legendre polynomial may be favoured with respect to the 
$ = uniform configuration, and thus mediate the tunneling events among the classical minima of the potential. 
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effect arising when the kinetic contribution to the action becomes negligible with respect to the 
potential one [86,87]. 

As figure shows, our results confirm that at the transition the average kinetic contribution is 
much smaller than the modulus of the potential. However, this behaviour persists (and is, in fact, 
enhanced) down to more strongly negative r-values in the uniform order regime, too — a fact 
that does not allow to interpret the existence of the striped phase as solely due to the potential 
dominance. 
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Figure 2: Average values of the kinetic term and modulus of the potential per degree of freedom, 
in the proximity of the transition between the disorder and striped phases. Data obtained at 
fixed A = |f. 



Rather, the results are compatible with the fact that the tunneling between the minima 
of the potential may be mediated by matrix configurations corresponding to non-spherically- 
symmetric distributions on the sphere. This can be easily justified from the analytical point of 
view (evaluating explicitly the action associated to such configurations), and is fully consistent 
with the numerical data (which confirm a non-negligible expectation value for the modes with 
I ^ when r is not very far from zero). 

Having presented the general features of the phase structure of the model, we now address the 
test of the theoretical predictions formulated in [33@: The model can be described by means of 
random matrix methods, and is characterised by well-defined properties of the eigenvalue sector. 
In particular, the phase transition is associated with a change in the topology of the support of 
the eigenvalue distribution. 



4 Note that the conventions in [33] differ with respect to eq. fT]) and to those used in [86] by a trivial rescaling 
of the $ matrix, and the coefficients of the quadratic and quartic terms in the potential are denoted as m and g, 
respectively. 
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We concentrated our attention onto matrices of size N = 15, 19, 21, 23, 31, 41 and 81. 

At intermediate and large values of g the data agree well with the theoretical prediction, and 
the theoretical prediction is indeed confirmed in the double-scaling limit. Figure shows the 
evolution of the eigenvalue density, from the one-cut to the two-cut phase. 

We should emphasise that the agreement observed is highly non-trivial, because the argu- 
ments underlying the derivation in [33] are expected to hold only approximately in two dimen- 
sions!! Although it would not be completely justified to perform a fit of the data to the curve, 
the agreement between the numerical results and the theoretical curve, which is a completely 
parameter-free prediction and does not involve any fitting procedure, is striking. 

On the contrary, the model exhibits significant quantitative deviations from the expected 
critical point at small g values. Although the qualitative pattern of the transition is confirmed, 
the location of the transition point is shifted towards more negative r values. In fact, this is 
not surprising: in this regime, the secondary maxima in the eigenvalue distribution are no longer 
negligible; the system could rather be described by a multi-trace model [110]. 

However, the general agreement of numerical results and theoretical predictions is also con- 
firmed by the scaling properties of the eigenvalue density support. The samples where the ob- 
served width of the distribution support differs from the expected one by an amount of the order 
of 5 - 10% or more correspond to strongly negative r values — which actually may even already 
lie in the "uniformly ordered" phase]! 

Table Q] shows results for the numerical critical values of rR 2 , for various matrix sizes. 

In the remaining part of this section, we shortly discuss the case when the classical potential 
has a unique minimum, and propose a method to detect the effect of the non-commutative 
anomaly. 

A perturbative study of the model in this regime was presented in [690. As it was discussed 
above, the non-commutative anomaly shows up as a (mild) non-local effect, distorting the energy- 
momentum relation on the fuzzy sphere by a finite amount. 

The one- loop effective action on the fuzzy sphere is [69]: 



while the non-commutative anomaly is given by the contribution involving h yAj ; h(x) is the 
harmonic number: h(x) = Y^t=i Ti with h(0) = 0, and A is a function of the Laplace operator, 

5 This is due to the weaker dominance of planar diagrams over the non-planar ones in D = 2. 

6 The critical line corresponding to the transition from the non-uniform to the uniform order phase was obtained 
numerically in [87], but a full theoretical description of this transition in the fuzzy sphere setting is still missing. 

7 In the notations of [69], the radius is set to the unit, the <E> matrix (whose size is denoted as TV + 1) is rescaled 
by a factor v2 with respect to eq. (JTJ, and the coefficients of the quadratic and quartic terms in the potential are 
denoted as fi 2 and Jj, respectively. 





where Sfj? is the square-mass renormalisation: 




(6) 
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(m') 4 / (4g) = 0.464 (m') 4 / (4g) = 1 .05 
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Figure 3: Evolution of the eigenvalue distribution from the one-cut (top-left plot) to the two-cut 
phase (bottom-right). The theoretical prediction for the transition point is: (m') 4 / (4</) = 1. 
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Table 1: Numerical and theoretical values for the quadratic coupling at the critical point. 



whose eigenvalue when acting on Yj m is I. In particular, the action of h ( A ) on the <5 matrix 
appearing in eq. (|2|) reads: 



N-l 



1=1 L \t=l / m=-l 



(7) 



In order to disentangle among the various terms contributing to the effective action, one can tune 
the parameter in such a way, that the main momentum-dependent contribution at order g comes 
from the non-commutative anomaly term only. Then, the spectrum of relative weights associated 
with the various spin channels gets distorted, as g is changed: larger probabilities are expected 
for higher I channels, when g is increased. 

Preliminary numerical tests (see figure H] as an example) confirm that the scalar component 
does not depend on g, while Hi(g), the power of the modes associated to higher I values, which 
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Figure 4: At large /i values, the channels associated to non-vanishing momenta are expected 
to be slightly enhanced by the non-commutative anomaly. The plot shows preliminary results 
from tests over an ensemble of matrices of size N + 1 = 7, at fi 2 = 200, in comparison with the 
theoretical prediction (dashed line). 

can be defined in terms of the coefficients appearing on the right-hand side of eq. ([2]) as: 



slowly increases with g, in rough agreement with the theoretical expectation. However, due to 
the theoretical approximations and technical difficulties involved in the observation of these fine 
effects, our data do not allow to make conclusive statements about this issue. This problem may 
be addressed more thoroughly in future work. 

4 Discussion 

The data in the previous section confirm the theoretical predictions for a scalar theory on the 
fuzzy sphere, as well as previous numerical results for the same model. The novel algorithm we 
used for the simulation proved to be very efficient, enabling us to obtain high-precision results for 
both the regimes that were investigated. The algorithm strongly reduces the correlation among 
subsequent configurations in a Markov chain, in a way which is compatible with the dynamics of 
the quantum system, and limiting undesired numerical artefacts to a minimum. 

As it concerns the case of a classical scalar potential with two degenerate minima, our data, 
obtained from simulations at a large number of points in the space of physical parameters of 
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the model, confirm and generalise the results obtained in similar numerical works [86,87]. We 
have discussed the role of non- uniform configurations on a fixed, finite radius sphere and their 
relevance to tunneling events connecting the classical vacua — a phenomenon which is not due 
to the non-commutative nature of the fuzzy regularisation, nor can it be directly interpreted as 
a signature of the UV/IR mixingH 

Yet, the relation to UV/IR mixing shows up, once one considers the double-scaling limit, 
and looks at the distribution for the matrix eigenvalues, which behave as a set of collective, 
intrinsically non-local degrees of freedom, and whose statistical properties can be worked out via 
random matrix methods. The UV/IR mixing manifests itself as the high-energy modes suppress 
the distribution of the low-energy ones. 

In our numerical study, we have successfully compared the observed eigenvalue distribution 
with the predicted pattern. The eigenvalues, rescaled through a factor predicted by the theory, fit 
very well to the [—1,1] range, and their density undergoes a transition from the "one-cut" to the 
"two-cut" regime. Although for the two-dimensional case one would only expect an approximate 
agreement, the data for p(ip) follow the theoretical curve very closely — except at very small 
^-values, where a multi-trace model [110] would probably provide a better description of the 
system. 

Next, we have also considered the m 2 > regime, and proposed an approximate method 
to observe the non-commutative anomaly, through the spectrum distortion at high momenta. 
Unfortunately, the parameter range in which unambiguous results can be obtained is severely 
limited, due to theoretical and numerical constraints. The results of a pilot study using this 
method look compatible with the expected effect, but for the moment they do not allow us to 
draw definite conclusions about the large- N limit; this aspect may be addressed again in future 
work. 

In conclusion, we can say that the non-perturbative results obtained in the present work 
confirm the current theoretical understanding of this simple non-commutative model, and the 
virtues and limits of the fuzzy approach as a regularisation scheme. 

On one hand, it is clear that the effects associated with non-commutativity are intrinsic to the 
fuzzy regularisation, and the presence of the anomaly discussed above threatens the possibility 
to use fuzzy spaces as a straightforward regularisation for QFT in ordinary (i.e. commutative) 
spaces. In this perspective, it would be particularly interesting to study in more detail the 
proposals that have been formulated [94, 95] to define an improved formulation of the action, 
yielding the correct QFT limit. 

On the other hand, the results discussed in this paper provide evidence that fuzzy spaces 
are indeed a well-suited regularisation scheme for theoreies directly defined in non-commutative 
spaces — e.g. Groenewold-Moyal Wg spaces — and offer the possibility for a practical and efficient 
implementation of numerical studies of these models. 

The success of this numerical study in the two-dimensional setting is very encouraging, and 
in the future it may be very interesting to address the D = 4 case. This generalisation would 

8 It may be instructive to point out that the situation for a finite-radius sphere is different with respect to the 
case of an infinite plane. 
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obviously be of great interest from the physical point of view, and — apart from a larger compu- 
tational effort — it is not expected to involve particularly difficult technical problems. In fact, 
the observation of the most interesting non-commutative effects may even turn out to be simpler 
than in D = 2 — due to the reasons discussed above. 
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A An overrelaxation algorithm for the fuzzy sphere 



In large parts of this work, a novel algorithm was used, that improves the efficiency of the 
numerical simulation. It allows to damp the autocorrelation between subsequent elements in the 
Markov chain of configurations. 

The basic idea is related to the overrelaxation technique in lattice gauge theory [101, 102]: 
The trial value in the update process of a given variable is chosen to be "as far as possible" from 
the original value, and such that the action of the system is unchanged^ In the lattice setting, 
this can be accomplished through a group reflection which can be worked out exactly for the 
SU(2) group, and in an efficient way for a generic SU(N) group [111,112]. This technique cannot 
be directly implemented in the present case, due to the fact that the $ variable takes values in a 
domain of different nature: the space of hermitian matrices of size N is non-compact, and, more 
important, a naive "reflection" of the $ matrix would not be effective for the purpose of reducing 
the autocorrelation time, since it would not allow to explore all of the physical orbits. 

The algorithm we built generalises the principia underlying the overrelaxation technique, 
adapting them to the present case, and it works as follows: Assume <l>o to be the initial matrix 
configuration, obtained with some ergodic procedure; let So = S($q) be the associated euclidean 
action. Let be a new, completely random (and, therefore, completely independent from <3? 
hermitian matrix in Mat at, with S* = S(<I>*) the corresponding value of the action. If S* > So 
a new hermitian matrix <J>i, such that: Si = S(<E>i) = So, can be built rescaling as: 



$i = (A.l) 



provided the following condition: 



|(S >0) V ^tr(**foj^,$*]]+ri2 2 $*) < -iJ^-^^tr^j j (A . 2) 

is true. If that is not the case, then is redefined (possibly iteratively) as: 

*. - ^ (A.3) 

until the condition in eq. (|A.2|) is satisfied. Note that this shift would drive <3?* closer and closer to 
$0) thus inducing a correlation between corresponding matrix entries in $i and $oj nevertheless, 
in general the eventual value obtained for <I>i may belong to a different physical orbit with respect 
to <J> . 

The algorithm is efficient under general conditions, including the cases in which S(<I>) is a 
function which varies strongly even for moderate changes in its argument, because the process 
driving towards $0 is exponentially fast, its implementation only involves trivial numerical 
operations, and terminates in a finite (and typically small) number of steps f^l 



9 This implies that the overrelaxation procedure is microcanonical; therefore the method is always combined 
with other canonical techniques, in order to ensure ergodicity of the whole update process. 
10 If S* < So, then is accepted as the new matrix configuration. 

11 This is easily proven using continuity and the fact that the trivial a — 1 solution exists for = <E>o- 
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500 1000 1500 2000 500 1000 1500 200O 500 1000 1500 2000 



Figure 5: Left: A typical Monte Carlo history of the action per degree of freedom and of tr(3>), 
obtained without overrelaxation; the data fluctuate close to one of the two minima of the potential, 
and tunneling events are very rare. Data obtained from matrices of size N = 27, in the two-cut 
regime. Center: Same as in the previous plot, but invoking the overrelaxation algorithm at Monte 
Carlo times r = 400, 800, 1200 and 1600. Right: Same as in the previous plot, but replacing the 
overrelaxation steps with new starts from unthermalised configurations. 

This algorithm proves superior to standard Metropolis (because it is not affected by the 
ergodicity problem) and to repeated new starts, because full thermalisation takes longer; figure [5] 
shows a comparison. 

The initial hermitian matrix <3?* can be chosen according to an arbitrary distribution; in order 
to achieve the best efficiency, we tested various possibilities, and in the production runs the matrix 
elements of were chosen according to a gaussian distribution centered around zero; the width 
of the gaussian is tuned according to an optimisation criterion. 

The number of ergodic updates between the overrelaxation steps is another tunable parameter 
of the algorithm. Typically, for the results presented here, one overrelaxation step was invoked 
every 50, 200, 500, 1000, or 2000 steps. 

References 

[1] S. Doplicher, K. Fredenhagen and J. E. Roberts, Commun. Math. Phys. 172 (1995) 187 
[arXiv:hep-th/0303037]. 

[2] G. Landi, arXiv:hep-th/9701078. 

[3] M. R. Douglas and N. A. Nekrasov, Rev. Mod. Phys. 73, 977 (2001) [arXiv:hep-th/0106048]. 

[4] R. J. Szabo, Class. Quant. Grav. 23 (2006) R199 [arXiv:hep-th/0606233]. 

[5] W. Pauli, "Scientific Correspondence," Vol. II, (Ed. Karl von Meyenn, Springer- Verlag, 1985), 
p. 15. 



14 



[6] W. Pauli, "Scientific Correspondence," Vol. Ill, (Ed. Karl von Meyenn, Springer- Verlag, 1993), 
p. 380. 

[7] H. S. Snyder, Phys. Rev. 71 (1947) 38. 

[8] C. N. Yang, Phys. Rev. 72 (1947) 874. 

[9] J. E. Moyal, Proc. Cambridge Phil. Soc. 45 (1949) 99. 

[10] D. Karabali, V. P. Nair and S. Randjbar-Daemi, arXiv:hep-th/0407007. 

[11] S. Bourouaine and A. Benslama, arXiv:hep-th/0610256. 

[12] A. Connes, M. R. Douglas and A. S. Schwarz, JHEP 9802, 003 (1998) [arXiv:hep- 
th/9711162]. 

[13] M. R. Douglas and C. M. Hull, JHEP 9802, 008 (1998) [arXiv:hep-th/9711165]. 

[14] A. Y. Alekseev, A. Recknagel and V. Schomerus, JHEP 9909 (1999) 023 [arXiv:hep- 
th/9908040]. 

[15] N. Seiberg and E. Witten, JHEP 9909 (1999) 032 [arXiv:hep-th/9908142]. 
[16] S. Kar, Nucl. Phys. B 577 (2000) 171 [arXiv:hep-th/9911251]. 
[17] S. Kar, Mod. Phys. Lett. A 18 (2003) 1053 [arXiv:hep-th/0006073]. 
[18] F. Lizzi, Int. J. Mod. Phys. A 22 (2007) 1317 [arXiv:hep-th/0610023]. 

[19] A. H. Chamseddine, G. Felder and J. Frohlich, Commun. Math. Phys. 155 (1993) 205 
[arXiv:hep-th/9209044]. 

[20] F. Lizzi, G. Mangano, G. Miele and M. Peloso, JHEP 0206 (2002) 049 [arXiv:hep- 
th/0203099]. 

[21] S. Doplicher, K. Fredenhagen and J. E. Roberts, Phys. Lett. B 331 (1994) 39. 

[22] S. Doplicher, arXiv:hep-th/0105251. 

[23] T. Filk, Phys. Lett. B 376, 53 (1996). 

[24] G. H. Chen and Y. S. Wu, arXiv:hep-th/0103020. 

[25] A. Jain and S. D. Joglekar, Int. J. Mod. Phys. A 19, 3409 (2004) [arXiv:hep-th/0307208]. 

[26] M. Chaichian, P. P. Kulish, K. Nishijima and A. Tureanu, Phys. Lett. B 604, 98 (2004) 
[arXiv:hep-th/0408069]. 

[27] R. J. Szabo, Phys. Rept. 378, 207 (2003) [arXiv:hep-th/0109162]. 



15 



[28] S. Minwalla, M. Van Raamsdonk and N. Seiberg, JHEP 0002, 020 (2000) [arXiv:hep- 
th/9912072]. 

[29] I. Hinchliffe, N. Kersting and Y. L. Ma, Int. J. Mod. Phys. A 19, 179 (2004) [arXiv:hep- 
ph/0205040]. 

[30] S. S. Gubser and S. L. Sondhi, Nucl. Phys. B 605, 395 (2001) [arXiv:hep-th/0006119]. 
[31] P. Castorina and D. Zappala, Phys. Rev. D 68 (2003) 065008 [arXiv:hep-th/0303030]. 
[32] E. Langmann and R. J. Szabo, Phys. Lett. B 533 (2002) 168 [arXiv:hep-th/0202039]. 
[33] H. Steinacker, JHEP 0503 (2005) 075 [arXiv:hep-th/0501174]. 

[34] H. Grosse and R. Wulkenhaar, Commun. Math. Phys. 256 (2005) 305 [arXiv:hep- 
th/0401128]. 

[35] V. Rivasseau, F. Vignes-Tourneret and R. Wulkenhaar, Commun. Math. Phys. 262 (2006) 
565 [arXiv:hep-th/0501036]. 

[36] T. R. Govindarajan, S. Kurkciioglu and M. Panero, Mod. Phys. Lett. A 21 (2006) 1851 
[arXiv:hep-th/0604061]. 

[37] H. Grosse and H. Steinacker, Nucl. Phys. B 746 (2006) 202 [arXiv:hep-th/05 12203]. 

[38] F. Lizzi, S. Vaidya and P. Vitale, Phys. Rev. D 73 (2006) 125020 [arXiv:hep-th/0601056]. 

[39] A. Stern, Nucl. Phys. B 745 (2006) 236 [arXiv:hep-th/0602062]. 

[40] H. Grosse and H. Steinacker, JHEP 0608 (2006) 008 [arXiv:hep-th/0603052]. 

[41] J. M. Gracia-Bondi'a, F. Lizzi, F. R. Ruiz and P. Vitale, Phys. Rev. D 74 (2006) 025014 
[Erratum-ibid. D 74 (2006) 029901] [arXiv:hep-th/0604206]. 

[42] H. Grosse and H. Steinacker, arXiv:hep-th/0607235. 

[43] M. Disertori and V. Rivasseau, arXiv:hep-th/0610224. 

[44] H. Grosse and M. Wohlgenannt, arXiv:hep-th/0703169. 

[45] J. Ambj0rn, Y. M. Makeenko, J. Nishimura and R. J. Szabo, JHEP 9911 (1999) 029 
[arXiv:hep-th/9911041]. 

[46] J. Ambj0rn, Y. M. Makeenko, J. Nishimura and R. J. Szabo, Phys. Lett. B 480 (2000) 399 
[arXiv:hep-th/0002158]. 

[47] J. Ambj0rn, Y. M. Makeenko, J. Nishimura and R. J. Szabo, JHEP 0005 (2000) 023 
[arXiv:hep-th/0004147]. 



16 



[48] T. Azuma, S. Bal, K. Nagao and J. Nishimura, JHEP 0405, 005 (2004) [arXiv:hep- 
th/0401038]. 

[49] W. Bietenholz, F. Hofheinz and J. Nishimura, JHEP 0209 (2002) 009 [arXiv:hep- 
th/0203151]. 

[50] W. Bietenholz, F. Hofheinz and J. Nishimura, Nucl. Phys. Proc. Suppl. 119 (2003) 941 
[arXiv:hep-lat/0209021]. 

[51] J. Ambj0rn and S. Catterall, Phys. Lett. B 549, 253 (2002) [arXiv:hep-lat/0209106]. 

[52] W. Bietenholz, F. Hofheinz and J. Nishimura, JHEP 0406 (2004) 042 [arXiv:hep- 
th/0404020]. 

[53] W. Bietenholz, A. Bigarini, F. Hofheinz, J. Nishimura, Y. Susaki and J. Volkholz, Fortsch. 
Phys. 53 (2005) 418 [arXiv:hep-th/0501147]. 

[54] W. Bietenholz, J. Nishimura, Y. Susaki and J. Volkholz, JHEP 0610 (2006) 042 [arXiv:hep- 
th/0608072]. 

[55] J. Madore, Class. Quant. Grav. 9 (1992) 69. 

[56] G. Alexanian, A. P. Balachandran, G. Immirzi and B. Ydri, J. Geom. Phys. 42, 28 (2002) 
[arXiv:hep-th/0103023]. 

[57] A. P. Balachandran, B. P. Dolan, J. H. Lee, X. Martin and D. O'Connor, J. Geom. Phys. 
43, 184 (2002) [arXiv:hep-th/0107099]. 

[58] A. B. Hammou, M. Lagraa and M. M. Sheikh- Jabbari, Coherent state induced star- 
product on R(lambda)**3 and the fuzzy sphere, Phys. Rev. D 66 (2002) 025025 [arXiv:hep- 
th/0110291]. 

[59] J. Medina and D. O'Connor, JHEP 0311, 051 (2003) [arXiv:hep-th/0212170]. 

[60] H. Grosse, C. Klimci'k and P. Presnajder, Commun. Math. Phys. 185 (1997) 155 [arXiv:hep- 
th/9507074]. 

[61] H. Grosse and G. Reiter, Jour. Geom. Phys. 28 (1998) 349 [arXiv:math-ph/9804013]. 

[62] A. P. Balachandran, S. Kurkciioglu and E. Rojas, JHEP 0207 (2002) 056 [arXiv:hep- 
th/0204170]. 

[63] S. Vaidya and B. Ydri, Nucl. Phys. B 671, 401 (2003) [arXiv:hep-th/0305201]. 
[64] B. P. Dolan and D. O'Connor, JHEP 0310 (2003) 060 [arXiv:hep-th/0306231]. 
[65] F. Lizzi, P. Vitale and A. Zampini, JHEP 0308 (2003) 057 [arXiv:hep-th/0306247]. 



17 



[66] F. Lizzi, P. Vitale and A. Zampini, JHEP 0509 (2005) 080 [arXiv:hep-th/0506008]. 

[67] A. P. Balachandran, S. Kurkciioglu and S. Vaidya, arXiv:hep-th/0511114. 

[68] M. M. Sheikh-Jabbari, Phys. Lett. B 642 (2006) 119 [arXiv:hep-th/0605110]. 

[69] C. S. Chu, J. Madore and H. Steinacker, JHEP 0108 (2001) 038 [arXiv:hep-th/0106205]. 

[70] S. Baez, A. P. Balachandran, B. Ydri and S. Vaidya, Commun. Math. Phys. 208 (2000) 787 
[arXiv:hep-th/9811169]. 

[71] H. Grosse, C. Klimcfk and P. Presnajder, Commun. Math. Phys. 178, 507 (1996) [arXiv:hep- 
th/9510083]. 

[72] A. P. Balachandran and S. Vaidya, Int. J. Mod. Phys. A 16, 17 (2001) [arXiv:hep- 
th/9910129]. 

[73] P. Presnajder, J. Math. Phys. 41, 2789 (2000) [arXiv:hep-th/9912050]. 

[74] U. Carow-Watamura and S. Watamura, Commun. Math. Phys. 212 (2000) 395 [arXiv:hep- 
th/9801195]. 

[75] H. Steinacker, Nucl. Phys. B 679 (2004) 66 [arXiv:hep-th/0307075]. 

[76] W. Behr, F. Meyer and H. Steinacker, JHEP 0507 (2005) 040 [arXiv:hep-th/0503041]. 

[77] V. P. Nair, Nucl. Phys. B 750 (2006) 321 [arXiv:hep-th/0605008]. 

[78] P. Valtancoli, J. Phys. A 39 (2006) 6069. 

[79] B. P. Dolan, I. Huet, S. Murray and D. O'Connor, arXiv:hep-th/0611209. 

[80] S. Murray and C. Samann, arXiv:hep-th/0611328. 

[81] C. Samann, arXiv:hep-th/0612173. 

[82] H. Aoki, S. Iso and T. Maeda, arXiv:hep-th/0610125. 

[83] R. Delgadillo-Blando and B. Ydri, JHEP 0703 (2007) 056 [arXiv:hep-th/0611177]. 
[84] H. Steinacker and R. J. Szabo, arXiv:hep-th/0701041. 

[85] H. Grosse, C. Klimcfk and P. Presnajder, Int. J. Theor. Phys. 35 (1996) 231 [arXiv:hep- 
th/9505175]. 

[86] X. Martin, JHEP 0404, 077 (2004) [arXiv:hep-th/0402230]. 

[87] F. Garcia Flores, D. O'Connor and X. Martin, PoS LAT2005 (2006) 262 [arXiv:hep- 
lat/0601012]. 



18 



[88] J. Medina, W. Bietenholz, F. Hofheinz and D. O'Connor, PoS LAT2005, 263 (2006) 
[arXiv:hep-lat/0509162]. 

[89] T. Azuma, S. Bal, K. Nagao and J. Nishimura, JHEP 0407, 066 (2004) [arXiv:hep- 
th/0405096]. 

[90] K. N. Anagnostopoulos, T. Azuma, K. Nagao and J. Nishimura, JHEP 0509, 046 (2005) 
[arXiv:hep-th/0506062]. 

[91] T. Azuma, S. Bal, K. Nagao and J. Nishimura, JHEP 0509, 047 (2005) [arXiv:hep- 
th/0506205]. 

[92] D. O'Connor and B. Ydri, JHEP 0611 (2006) 016 [arXiv:hep-lat/0606013]. 
[93] M. Panero, SIGMA 2 (2006), 081 [arXiv:hep-th/0609205]. 

[94] B. P. Dolan, D. O'Connor and P. Presnajder, JHEP 0203 (2002) 013 [arXiv:hep-th/0109084]. 
[95] B. P. Dolan and D. O'Connor, private communication. 

[96] M. L. Mehta, "Random Matrices," Pure and Applied Mathematics, Volume 142, 3rd edition 
(Academic Press, San Diego, 2004). 

[97] M. A. Stephanov, J. J. M. Verbaarschot and T. Wettig, "Random matrices," originally 
published in the Wiley Encyclopedia of Electrical and Electronics Engineering, Supplement 1 
(2001) [arXiv:hep-ph/0509286]. 

[98] E. Brezin, C. Itzykson, G. Parisi and J. B. Zuber, Commun. Math. Phys. 59 (1978) 35. 

[99] Y. Shimamune, Phys. Lett. B 108 (1982) 407. 

[100] P. Bleher and A. Its, arXiv:math-ph/0201003. 

[101] S. L. Adler, Phys. Rev. D 23 (1981) 2901. 

[102] F. R. Brown and T. J. Woch, Phys. Rev. Lett. 58 (1987) 2394. 

[103] A. D. Kennedy and B. J. Pendleton, Phys. Lett. B 156 (1985) 393. 

[104] K. Fabricius and O. Haan, Phys. Lett. B 143 (1984) 459. 

[105] M. Liischer and P. Weisz, JHEP 0109 (2001) 010 [arXiv:hep-lat/0108014]. 

[106] M. Zach, M. Faber and P. Skala, Phys. Rev. D 57 (1998) 123 [arXiv:hep-lat/9705019]. 

[107] M. Panero, JHEP 0505 (2005) 066 [arXiv:hep-lat/0503024]. 

[108] N. Madras and A. D. Sokal, J. Statist. Phys. 50 (1988) 109. 



19 



[109] J. Shao and D. Tu, "The Jackknife and Bootstrap," (Springer, New York, 1995). 

[110] D. O' Connor and C. Samann, in preparation. 

[Ill] R. Petronzio and E. Vicari, Phys. Lett. B 248 (1990) 159. 

[112] P. de Forcrand and O. Jahn, arXiv:hep-lat/0503041. 



20 



